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Abstract 

We study a quantum Hamiltonian that models an two-dimensional array 
of Josephson junctions with short range Josephson couplings, (given by the 
Josephson energy Ej) and charging energiey Eq due to the small capacitance 
of the junctions. We include the effects from both the self-Co and the junction- 
C\ capacitances in the presence of external magnetic flux / = ^/^o as well 
as uniform background of charges q x . We derive an effective quantum non- 
linear (T-model for the array Hamiltonian which enables us a non mean-field 
treatment of the zero-temperature phase transition scenario. We calculate the 
ground-state phase diagram, analytically deriving Ej nt (Ec,q x , f) for several 
rational fluxes / = 0,^,|,| and ^ that improves upon previous theoretical 
treatments based on mean-field approximations. 
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I. INTRODUCTION 



Advances in the photolithographic micro-fabrication techniques have experimentally al- 
lowed to manufacture Josephson-j unction arrays (JJA) with tailor specific properties. Be- 
cause the junction parameters can be accurately controlled JJA offer a unique opportunity 
to test different quantum mechanical models. 1 In JJA's the two main energy scales are set by 
the Josephson coupling, Ej, between superconducting islands due to Cooper pair tunneling, 
and to the electrostatic energy E c (proportional to the inverse of capacitance matrix setup 
by island self- C and junction- C\ capacitances) due to local deviations from charge neutral- 
ity. When Ej is much larger then E c the phases of the superconducting order parameters on 
individual islands are well defined. In this regime the semi-classical fluctuations determine 
the properties of the JJA. In the opposite limit, ie. for Ec ^> Ej, the superconducting 
phases are dominated by strong quantum fluctuations since the Coulomb blockade pins the 
charge carriers to the islands. Several theoretical 2-4 and experimental studies 5 have con- 
sidered the competition between the Ej dominated phase and the E c dominated charging 
energy regions in periodic JJA. It has been established that for sufficiently large charging 
energy the quantum phase fluctuations lead a complete suppression of long-range phase 
coherence even at zero temperature. This type of quantum phase transition has attracted 
significant interest in recent years (for a review see Ref. 6). 

Because the junction parameters can be controlled accurately JJA offer a unique system 
where one can test the nature of quantum phase transitions and critical phenomena, in par- 
ticular the superconductor-insulator (SI) phase transition induced by quantum fluctuations. 
The experimental systems can be modeled by a quantum generalization of the classical XY 
model with quantum phase of the superconducting order parameter on each island canoni- 
cally conjugated to the excess of Cooper pair number. The quantum parameter a = E c /Ej 
determines the relevance of the quantum fluctuations. Bulk of theoretical studies have been 
carried out within the self-capacitance model (Co 7^ 0, C\ = 0), usually using different kinds 
of mean-field-theory or self-consistent harmonic approximations (for a comprehensive ac- 
count see, Ref. 1). In the mutual-capacitance dominated limit [C\ 7^ 0, Co = 0) the lattice 
model is equivalent to a quantum Coulomb gas model, 7 with critical properties that are not 
fully understood at present. Further complication arise when both C and C\ are non-zero. 
The charging energy can also be tuned by an applied gate voltage V g which control the 
energy required to change the number of Cooper pairs on the islands. A uniform V g can be 
applied by means of gate capacitance C g . 8 ' 9 Assuming C g <C Co the applied voltage induces 
a charge C g V g on each grain which can be conveniently account for when one defines the a 
dimensionless charge frustration parameter q x = C g V g /2e. The effects of charge frustration 
is closely related to the so-called Bose-Hubbard model 10 which describes the strongly in- 
teracting bosons under the competing kinetic energy (related to Josephson energy) and the 
potential energy (i.e charging energy in the JJA) effects. The on-site (chemical) potential 
corresponds in this analogy to the charge frustration parameter q x in JJA. Frustration in 
quantum JJA can be also introduced by applying magnetic field B. The presence of B 
induces vortices in the system described by the magnetic frustration parameter / (= $/<E»o> 
where <£> is the magnetic flux piercing a 2D lattice plaquette and $ = h/2e stands for the 
elementary flux quantum). Of special interest are cases when f = p/q with p and q being 
the rational numbers. Particularly prominent is the case of / = 1/2 because of the interplay 
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between continuous U(l) symmetry of the superconducting parameter and discrete Ising 
Z 2 symmetry associated with antiferromagnetic arrangements of the plaquette chiralities 
describing directions of the circulating currents in each plaquette. 11 

Usually, mean-field calculations or semiclassical approaches are not expected to be re- 
liable at T = and be capable to handle spatial and quantum fluctuation effects properly 
especially in two-dimensions. Notably, there are virtually no works beyond mean-field level 
on 2D JJA in the where the physics is dominated by the by the combination of charge q x , 
magnetic frustration / and the charging energy related to the self-Co and junction-Ci capac- 
itances. Recent studies based on perturbative expansions combine either JJA with charge 
frustration and and charging energy effects 12 or magnetic frustration and effects of array 
capacitances. 13 Therefore, the problem we would like to address is then: What is the effect 
of having a competition between the magnetic, charge frustration and quantum effects on 
the ground state ordering in the Josephson-j unction network on a square two-dimensional 
lattice? The purpose of this work is therefore to study these quantum fluctuation effects in 
an analytical way to understand the ground state phase diagram of the system as a function 
of various control parameters like q x and / as well as the ratio of the junction to the self- 
capacitances Ci/Cq. To analyze the 2D JJA beyond the mean-field theory we employ the 
path-integral formulation explicitly tailored for the macroscopic JJA Hamiltonian. The ef- 
fective action formalism allows then for an explicit implementation of the magnetic field and 
offset charges effects with the nearest-neighbor mutual and self-capacitances included. Fur- 
thermore, we adopt an approach based on the quantum phase fluctuation algebra, to map 
the microscopic JJA model Hamiltonian onto an effectively constrained system - a solv- 
able quantum spherical model in two- dimensions that captures both quantum and spatial 
fluctuation effects beyond the mean field level. 

There are four sections in this paper: Section II is devoted for the introduction of the 
microscopic Hamiltonian of the 2D JJA with a set of parameters describing with interactions 
and set of control parameters. In Sec. Ill the system is transformed into quantum action of 
the XY model under the magnetic field using the "imaginary-time" Matsubara approach. 
Furthermore, Sec. IV presents the resulting phase diagrams and finally, Sec. V summarizes 
the results and presents some discussions. In Appendix we gave a derivation of closed- 
form analytical formulae for the density of states of a two-dimensional square lattice with 
magnetic flux for a number of rational values p/q which might be also of general interest. 

II. QUANTUM PHASE HAMILTONIAN 

A Josephson junction array can be modeled by a periodic lattice of superconducting is- 
lands separated by insulating barriers. Each island becomes independently superconducting 
about the bulk transition temperature T c0 and it is characterized by an order parameter 
^(rj) = \^j (v i )\e 1 '^ , ' i \ where rj is a two-dimensional vector denoting the position of each 
island. The magnitude of the order parameter, |-?/>o( r i)|? is non-fluctuating when the temper- 
ature is lowered further and the onset of long range phase coherence due to the tunneling of 
Cooper pairs between the islands is responsible for the zero resistance drop in the arrays. The 
competition between the Josephson tunneling and the charging (Coulomb) energy, without 
dissipation, can be modeled by the Hamiltonian 
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TC— He + 7~ij, 

^ r 

Hj= Y, J(|ri-r 2 |)[l-cos(0 ri -0 r2 )]. (1) 

( r l. r 2> 

Here, Q r = (2e/i)d/d<f) r is the charge operator while r represents the superconducting phase 
operator of the grain at the site r; J(|ri — r 2 |) is the site-dependent Josephson coupling 
and [C _1 ] rr / is the inverse capacitance matrix. The first part of the action (8) than defines 
the electrostatic energy with C rr ' being a geometrical property of the array. This matrix is 
normally approximated (both theoretically and in experimental interpretation) as a diagonal 
(so-called self-capacitance Co) and mutual one C\ between nearest neighbors. In the case 
of an square 2D JJA of interest here we write the capacitance matrix as: 

C Tr > = (C + 4Ci) S rr/ - Ci Y <^r,r'+d (2) 

d 

with the vector d running over the nearest neighboring (n.n) islands. 

A. Effect of applied magnetic field 

A perpendicular magnetic field B given bt the vector potential A enters the Hamiltonian 
(1) through Peierls phase factor according to 

J(|r, - r,-|) - J B {\v t - r,|) = Jflr, - r,|) exp £ J A • dl) . (3) 

Thus, the phase shift on each junction is determined by the vector potential A and in a 
typical experimental situation it can be entirely ascribed to the external field B. We assume 
throughout this paper that the model (1) is defined on a square lattice with lattice spacing 
a. From Eq.(3) it follows that the properties if the array will be periodic with a period 
corresponding to the one flux quantum $ — hc/2e per plaquette. Of special interest are 
the values of the magnetic field which corresponds to the rational values of / = <&/<&o, 
e.g. / = 1/2, 1/4, 1/4, . . . where $ = Ba 2 . Since all properties of the Hamiltonian (1) are 
invariant under / — > / + 1 and also under / — > — / it is sufficient to consider / in the range 
< / < 1/2. 

B. Offset charges 

Offset charges are an important ingredient in the experimental array samples made of 
ultrasmall junctions. Several authors have shown that static background charges can have a 
pronounced effect on the SI transition at zero temperature. 14,15 Including the offset charges, 
q x , in the charging energy of Eq. (8) gives 

~\ E[ C_1 ]rr' (4 - q x ) (&, - q x ) . (4) 

r 

Thus, offset charges, or an external gate voltage applied between the array and the substrate, 
behave like a chemical potential for injection of Cooper pairs into the array. 
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III. JJA "PHASE ONLY" ACTION 



It is useful to derive a field-theoretic representation of the partition function for Eq (1). 
A convenient procedure is to introduce a path-integral representation in a basis diagonal in 
<f>j. Given the Hamiltonian Tl, one can simply perform a Legendre transformation to obtain 
the corresponding Euclidean Lagrangian in the Matsubara "imaginary time" r formulation 
(0<T<l/k B T = (3): 



(5) 



The partition function Z = Tre n / kBT of the system is then given by 



n / n 



DQr 

2tt 



■ J^drC[Q,4>] 



(6) 



Performing the integration over the charge variables the partition function is expressed as 



-/ 



(7) 



Here the functional integral is evaluated over the phases restricted to the compact interval 
[0, 2n] and with an effective "phase-only" action (h = 1) 



r,r' 



' d<j>r \ (dcjh 
dr { dr 



Sj[<l>]= fdr J2 *(|ri-r 2 |){l-cos[0 r ,(T)-0„(T)]}. 



(8) 



(ri,r 2 > 



Since the values of the phases 0, in Eq.(8) which differ by 2ir are equivalent, the path 
integral in Eq.(7) can be written in terms of non-compact phase variables 9 r (r), defined on 
the unrestricted interval [— oo, +oo], and by a set of winding numbers {n r } = 0, ±1, ±2, . . ., 
which are integers running from — oo to +oo (and physically reflects the discreteness of the 
charge 14 ), so that r (r) = r (O) + 2im r r/f3 + r (r). 



A. Effective non linear a model 

Bulk of existing analytical works on quantum JJA have employed variety of mean-field- 
like approximations 1 ' 16-18 which, unfortunately, are not reliable for treatment spatial and 
temporal quantum phase fluctuations in a low dimensional systems. To study the JJA 
model it appears at first natural to use a description in terms of an effective Ginsburg- 
Landau functional derived from the microscopic model of Eq. (8). Several studies of JJA 
have followed this route, also known as the coarse grained approach first developed by 
Doniach. 19 The essence of this method is to introduce a complex field order parameter ip T 
(or equivalently a two component real field) whose expectation value is proportional to 
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(exp(i0 r )). The non-zero value of this quantity describes the "phase-locking" or long-range 
phase ordering in the model. Unfortunately, the system is governed by the Ginzburg-Landau 
functional only as long as the order parameter is small. This is a serious shortcoming of the 
coarse grained approach since the method is restricted to the region of parameters in the 
vicinity of the critical point. It is therefore desirable to extend the coarse-grained method 
in a way which enables one a self-consistent description of the phase transitions in JJA 
in an extended parameter range. Following Ref. 20 we note that the model (8) encodes 
phase fluctuation algebra given by the Euclidean group E 2 constituted by the commutation 
relation between particle L r and phase (ladder) operators P r = e l<i>v : [L r ,P r r] = —PS rr i, 
[L T , P r t] = P^5 rr r and [P r , P r t\ = with the conserved quantity (invariant of the E 2 algebra) 



PrPl = 1 



(9) 



Thus, the proper theoretical treatment of quantum JJA must maintain the constraint (9). 
To proceed we write the partition function Z = J [Hi T>(pi} e~ s for the model (8) in terms 
path integral representation 21 by introducing the auxiliary complex fields ipi(r) which replace 
the original ladder operators Pi. Furthermore, relaxing the original "rigid" constraint (9) 
and imposing the weaker (spherical) condition J2iPrP} = N allows us to formulate the 
problem in terms of the (exactly) soluble quantum spherical (QS) model (see, Ref. 22). 
This can be conveniently done using the Fadeev-Popov method with Dirac delta-functional 
which facilitate both the change of integration variables and superimposition of the spherical 
constraint: 



S (El^(r)| 2 - Nj e 



■SjW 



x 



X 



e- 5 ^Ip[ReW(T) - W(r)] 



x 6\LmMT)-S»(<f>(T))]. 



(10) 



A convenient way to enforce the spherical constraint is the functional analog of the 
S— function representation 8(x) = j^(d\/2n)e lXx which introduces the Lagrange multi- 
plier A(r) thus adding an additional quadratic term (in -^—fields) to the action (8). The 
evaluation of the effective action in terms of the ip to the second order in ipi{r) gives the 
partition function of the quantum-spherical model Z = Zq S : 



Z QS - 



2tH 



-S QS [ip,\] 



(11) 



where the action of the effective non-linear cr-model reads 

5 QS [^,A]= f drdr'l £ [(J(| ri _ r2 | + A(r)5 ri , r2 )5(r-r') 

° l(ri,r 2 ) 

+ r 02 + (r 1 r; r 2 r')] <(r)^ r2 (r') - NX(r)5(r - r')} , 



(12) 



Furthermore, r o2 + (r 1 r; r 2 r') = [W / 02 1 ] + (r 1 T;r 2 T l ) is the two-point phase vertex correlator 
and 
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W 02 + (r 1 r;r 2 T' 



Z ° {n r } r 



2tt 



d6(0) 



S(0)+27rn r 
9(0) 



^We'^W-^^^- 50 ^, (13) 



with 



p2tt f e(0)+2wn r „ r/il 

{nr} r -A) ^(0) 



(14) 



where Z is the statistical sum of the "non-interacting" system described by the action 
Sc[0]. The phase-phase correlation function becomes: 14 



W 02 (ue) 



8E, 



c 



E 



exp 



-4:pE c (q - q x 



Z V (4^c) 2 - [8E c {q - q x ) - iu e 



2 J 



(15) 



where Z = J2 q ex P [~ ^(3Ec{q — q x ) 2 }, and the summation is performed over all integer- 
valued charge states q — 0, ±1, ±2, . . ., which makes the function W^i^e) periodic. At low 
temperatures the sum over q in Eq.(15) is dominated by the charge q which makes the 
exponent in the numerator of Eq.(15) smallest (for T = this value is q — 0). 

In the thermodynamic limit, N — > oo, we can calculate the functional integral in Eq.(ll) 
by the steepest-descent method. To proceed, we introduce propagators associated with the 
order parameter field defined by 

G(r 1 r;r 2 r') = «(r)V r2 (r')) QS . (16) 

The condition that the integrand in Eq.(ll) has a saddle point A(r) = A is that 



(17) 



which becomes an implicit equation for X . The QS ensemble average is now defined by 



(...) 



QS 



f[n T T>rl> T VP T 



. e 



I [Ilr VlprVlpf] e- 5 Qs[V>,Ao] 



(18) 



To proceed, it is desirable to introduce the density of states for a 2D lattice with the magnetic 
flux / = p/q in the form 



Pp/q( E ) = ^T, 5 



E — 



Ej 



(19) 



with Je(k) the Fourier transform of the Josephson interactions Js(ji"i — r 2|) (cf. Eq.(3)) in 
a magnetic field. The problem of computing of p p / q (E) reduces effectively to the solution of 
the Harper's equation 23 relevant e.g. for tight-binding electrons on a two-dimensional lattice 
with a a uniform magnetic flux per unit plaquette. 24 In Appendix we give an analytical 
derivation of p p / q (E) in a closed form for several rational values of p/q. A Fourier transform 
of Eq. (12) in momentum and frequency space enables one to write the spherical constraint 
(17) with the help of the density of sates (19) in the form 
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1=4 r>E — — — i' (2°) 

where W / o2(^) is the frequency transformed phase-phase correlator of Eq. (13) and u>e = 
2ir£/(3 (£ = 0, ±1, ±2, . . .) the (Bose) Matsubara frequencies. As usual in a spherical model 
the critical behavior and the phase transition boundary depends crucially on the spectrum 
given by the density of states p p / q (E) and is determined by the denominator of the summand 
in the spherical constraint equation of Eq. (17). Specifically, when [1/G(k = 0, uji — 0)] = 0, 
where G(k, u>e) is the Fourier transformed order parameter correlation of Eq. (16), the system 
displays a critical point at 

A - e™ f a *Ej + [W 02 (cu e = 0)]- 1 = 0, (21) 

where e™?* is the maximum value of the spectrum described by the density of states (19). 
This fixes the saddle point value A: with the onset of the phase transition saddle point value 
of the Lagrange multiplier A "sticks" to that value at criticality (A = A ) and stays constant 
in the whole low temperature phase. 



IV. PHASE DIAGRAMS 

By substituting the value of A from Eq.(21) into (20), after performing the summation 
over Matsubara frequencies, one obtains the T — > limit the result 



+°° E, 



-oo 



C 



2(A + 2E C - £Ej) 



sign I 4 <7*-|r + 



y/2{\ + 2E C - 


£Ej)E c \ 


Ej j 


y/2(X Q + 2E C 


- ^Ej)E c 



sign j , . 



, (22) 



where V denotes the principal value of the integral. The charging energy for the two- 
dimensional square lattice may be explicitly written as 

E c = l -e 2 [C-\ r 

t C l rr ~ n^oo N ^ C + 4Ci - 2Ci(cos k x + cos k y ) 

2 (23) 



7r(C + 4Ci) VC0 + 4C! 

where K(x) is the elliptic integral of the first kind 25 (see also Eq.(A5) in Appendix). Fur- 
thermore, we introduce the charging energy parameters Eq = ^6 2 Cq 1 and E\ = ^e 2 Ci 1 
related to the island and junction capacitances. We fist examine the ground state phase 
diagram for a number of rational values of the flux piercing the array without charge frus- 
tration (q x = 0). We give in Fig.l a graph showing the ground state phase diagram as a 
function of island-to-junction charging energy for various rational flux numbers. 
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The charge frustration q x may serve as a control parameter distinguishing between or- 
dered and non-ordered states. The energy difference for two charge states in each grain of 
the JJA with n and n + 1 extra electrons may be reduced by changing q x . As a consequence 
the effects of a finite charging energy E c are weakened and the superconducting region in the 
phase diagram turns out to be enlarged. The dependence of Josephson to the total charging 
energy Eq.(23) Ej/Eq is shown in Fig. 2 for a number of rational fluxes. The behavior shown 
in in Fig. 2 is readily understood by realizing that the Hamiltonian (1) is periodic in q x with 
period unity, this implies that the phase diagram should be similarly periodic, repeating at 
each integer number. 14 



V. SUMMARY 



We have studied the ground phase diagram in quantum two-dimensional Josephson 
junction arrays, using the spherical model approximation with exactly evaluated density 
of states for two-dimensional lattice with rational magnetic flux / = p/q for a number of 
values / = \jq. One of the most important issues in these systems is the role of the charging 
energy on the phase coherence transition. The detailed phase boundaries crucially depend 
on the ratio of the junction-to-self- capacitances, Ci/C as well as on the magnetic and 
charge frustration. In arrays which are in the superconducting state at / = but with the 
value Ec/Ej close to the critical value, a magnetic filed can be used to drive the array into 
the insulating state (see, Fig.l). It is interesting to compare the calculated values of Ec/ Ej 
for / = and / = 1/2 with the experiments: 5 the measured critical value E c /Ej close to 
T = is about a factor 0.7 lower than the zero field value. This is to be compared with the 
value 0.88 found for T = in the present paper. Note that the critical value of Eq/Ej is 
greater than for others commensurate values f — 1/q which suppress the superconducting 
state substantially. We found that the effect of this supression is considerably larger than 
that which comes from the mean-field analysis of JJA with magnetic flux (for comparision 
see, e.g. Ref. 14). Closing we note that a reliable comparison with experiments presumably 
will still require incorporating into the model several other ingredients as, e.g, disorder and 
dissipation. These are topics to be considered in the future. 



APPENDIX A: DENSITY OF STATES FOR 2D SQUARE LATTICE IN A 

RATIONAL FLUX CASE 



If one uses the Landau gauge, A = 5(0, x, 0) than the dispersion e(k) for a square lattice 
of spacing a with flux = _Ba 2 /$o = p/q with is given by: 24 



det 



/ Mi -e ikxa 

_ e -ik x a gik x a 

_ e -ik x a M 



_ e -ik x a \ 




— e 



ik x a 



o ••. ••. -< 

J,k T a n 



%krr CL 



M„_i -e 



Xkrr CL 



q . 

D ik x d 



(Al) 



9 



where 



M„ 



-2 cos (k y a + 27r0n) — e (k) . 



(A2) 



Equation (A2) is known as Harper's equation and has been studied extensively. If integers 
p and q are chosen to represent the flux (with no common factor in p and q), then the 
dependence on the wave vector always appears through the generalized structure factor 
7 n = cos(nk x a)+cos(nk y a). The density of states pt(E) given by the Eq.(19) can be obtained 
by computing energy bands e(k) from the eigenvalue equation (Al). The calculation of the 
exact density of states is straightforward, although for large values of q may only be done 
numerically. 24 However, for a number of q values of interest it can be done analytically with 
a closed-form expression for pz(E) as end result. Below we list these cases. 



1. f = 

In the case of zero magnetic field the density of states for the square two-dimensional 
lattice reads simply 



p(E) 



77 /a (fk 



5[£-7i(k)]. 



L/a (27r) 2 

Performing the integration over the wave-vectors we obtain (see, Fig. 3) 



(A3) 



P2D (E) = — K 



7T Z 





(A4) 



where 



K(x) 



t/2 d 



/o yjl - x 1 sin 2 
is the elliptic integral of the first kind 25 and 

' 1 for x > 



(A5) 



e(x) 



for x < 



(A6) 



the unit step function. 



2- / = * 



2 



For the square lattice with the half-flux quantum per plaquette the Harper's equation 
(A2) reads 

e 2 (k) - 2 72 - 4 = (A7) 
so that the energy dispersion has two branches: 

10 



e(k) 



-^/4 + 2(cos2fci + cos2fc 2 ) if E e [-2^2; o) 
+ ^4 + 2 (cos2A;i + cos2£;2) if -E e |0; 2y / 2 _ J 



(A8) 



Integrating over the wave vectors belonging to the Brillouin zone we obtain from Eq.(A3) 
with the help of the formula (A4) 



2 ' \ 2 
Fig. 4 gives the plot of the density (A9). 



(A9) 



3 - f = \ 

For p/q = 1/3 the Harper's equation is now of cubic form 

e 3 (k) - 6e(k) + 2 73 = 
and its roots give the dispersion of the energy 

l + iV3 

r 

3 



(A10) 



73 + V- 8 + 73 



+i(i + y/3) (-13 + ^^8 + 



if E e (-1 - v^; -2 



"73 



(All) 



2 + (-73 + \/- 8 + 



7l 



-\{l + iV3) (-73 + V / - 8 + 73 2 ) 3 if £e(l-V3; 

if E E (2; 1 + v 7 ^ 



■1 + V3 



-73 + V -8 + 73 



The density of states then becomes (see, Fig. 5) 



Pi (#) = ^P2D 



X 

+ 
+ 
+ 



E (E 2 - 6) 



2 


-1 


13 -2 





+ 2 

a: i +2 



£3 - 2 



-1 



(£ 2 - 8) (£ 2 - 2) £3- 

(E + 1 + v^) - 6 (E + 2) 
6 (E - 1 + V3) - (£)] + 

e (£?) - e (e - i + V3)] + 
e (£ - 2) - e (e - i - 



X 



+ 



where 



x 



(E) = ^(e 3 -QE + ^E 2 (E 2 - 6) 2 - 32^ . 



(A12) 

(A13) 
(A14) 
(A15) 
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The matrix formula (A2) produces now a four order polynomial equation for the disper- 
sion e(k). 

e 4 (k)-8e 2 (k)- 2 74 + 4 = (A16) 
Because the Eq.Q is in fact bi-quadratic with respect to e(k) its roots can be readily found: 

' -^4 + v/VSTtI if E e -2y/2^ 



J(k) = 



+JA - V2^WTYa if Ee (0;V4-2v^ 



(A17) 



+ ^4 + V^VQTYi H Ee (2V2; ^4 + 2^2 



As a result the density of states becomes 



£ 2 -4| {E 4 -8E 2 + 4\ 

r P2D z X 



x j ^4 + |£ 2 - 4| [6 (8 - E 2 ) - (4 + 2^2 - E 2 ) 



+ 



+ yj4- \E 2 -4\& (4 - 2V2 - E 2 
In Fig. 6 we have plotted the density of states given by Eq.(Al9). 

5- f = \ 

The Harper' s equation (A2) now is given by 

e 6 (k) - 12e 4 (k) + 24e 2 (k) - 2 7e - 4 = 



(A18) 



(A19) 



(A20) 



which is a bi-cubic equation for the dispersion parameter e(k). Although the solutions 
of Eq.(A20) can be found analytically the result is quite involved so that we refrain form 
reproducing it here. However, the resulting density of states (see, Fig. 7) can be put in a 
closed-form expression as follows: 



pi (E) = 2 sp 2D 



E 6 - 12E 4 + 24E 2 - 4 



x 



x < I — 

,x 



2\s y/x* + 4 ' ( 2x )^ + 8 - 2 ^ 



{2x 2 Y - 16 



(5 + V21 - E 2 ) - (6 + 2V3 - E 2 ) 



+ 



^4-2§ - (2x^j 


-1 


23 (l - 


xi + 16 




23 




-i>/3) 
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FIGURES 



FIG. 1. Zero temperature phase diagram for two-dimensional Josephson junction array on a 
square lattice in magnetic field as a function of charging energies related to self-capacitance (Eq) 
and mutual capacitance (E\), (normalized by the Josephson coupling energy Ej) for different 
values of the flux ratio /. The inset shows closeup of the plot for / = and / = 1/2. The 
superconducting (insulting) phase is located above (below) each curve. 

FIG. 2. Ground state phase diagram of the square J J A as a function of the total charging 
energy (Ec, normalized by the Josephson coupling Ej) and the charge frustration parameter q x . 
The inset shows closeup of the plot for / = and / = 1/2. The insulating phase is located within 
each lobe. 

FIG. 3. Density of states p{E) for two-dimensional square lattice. 

FIG. 4. Density of states for two-dimensional square lattice with one-half of the magnetic flux 
quantum per plaquette. 

FIG. 5. Density of states for 2D square lattice in the magnetic field for the flux ratio / = 1/3. 

FIG. 6. Density of states for 2D square lattice in the magnetic field for the flux ratio / = 1/4. 

FIG. 7. p(E) for 2D square lattice in the magnetic field for the flux ratio / = 1/6 (depicted for 
positive argument E only, to reveal details of the band structure. 
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